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ABSTRACT 

An extreme-mass-ratio burst (EMRB) is a gravitational wave signal emitted when a 
compact object passes through periapsis on a highly eccentric orbit about a much more 
massive object, in our case a stellar mass object about a 1O 6 M black hole. EMRBs are 
a relatively unexplored means of probing the spacetime of massive black holes (MBHs) . 
We conduct an investigation of the properties of EMRBs and how they could allow 
us to constrain the parameters, such as spin, of the Galaxy's MBH. We find that if 
an EMRB event occurs in the Galaxy, it should be detectable for periapse distances 
r p < 65r g for a \i — 1OM orbiting object, where r g = GM./c 2 is the gravitational 
radius. The signal-to-noise ratio scales as log(p) ~ — 2.71og(r p /r g ) -f log(/i/M ) +4.9. 
For periapses r p < 10r g , EMRBs can be informative, and provide good constraints on 
both the MBH's mass and spin. Closer orbits provide better constraints, with the best 
giving accuracies of better than one part in 10 4 for both the mass and spin parameter. 

Key words: black hole physics - Galaxy: centre - gravitational waves - methods: 
data analysis. 



1 BACKGROUND AND INTRODUCTION 



related to the BH's angular momentum J by 



Many, if not all, galactic nuclei have harboured a massive 
black hole (MBH) during their evolution ( |Lynden-Bell &; 
Rees 1971; Rees 1984). Observations show there exist well- 
defined correlations between the MBHs' masses and the 
properties of their host galaxies, such as bulge luminosity, 
mass, velocity dispersion and light concentration (e.g. 



Ko- 



rmendy & Richstone 1995; Magorrian et al. 1998; Gra 
et al.|2001||Tremaine et al.|2002||Graham et al.|2011|). These 



suggest coeval evolution of the MBH and galaxy ([Peng 2007 
[Jahnke &; Macc io 2011), possibly with feedback mechanisms 
coupling the two ( |Haiman &; Q uataert 2004| |Volonteri &; ' 
Natarajan 2009}. The MBH and the surrounding spheroidal 
component share a common history, such that the growth 
of one can inform us about the growth of the other. 

The best opportunity to study MBHs comes from the 
compact object in our own galactic centre (GC), which is 
coincident with Sagittarius A* (Sgr A*). Through careful 
monitoring of stars orbiting the GC, this has been identified 
as an MBH of mass M. 4.31 x 1O 6 M at a distance of 
only Rq = 8.33 kpc ( |Gillessen et al.|2009) . 

According to the no-hair theorem, the MBH should be 
described completely by just its mass M. and spin a, since 
we expect the charge of an astrophysical black hole to be 



negligible ( Chandrasekhar 1998). The spin parameter a is 



J M.ac; 

it is often convenient to use the dimensionless spin 
_ cJ 



(1) 



(2) 



As we have a good estimate of the mass, to gain a complete 
description of the MBH we have only to measure its spin; 
this shall give us insight into its history and role in the 
evolution of the Galaxy. 

The spin of an MBH is determined by several competing 
processes. An MBH accumulates mass and angular momen- 
tum through accretion (Volonteri 2010). Accretion from a 



gaseous disc shall spin up the MBH, potentially leading to 
high spin values ( Volonteri et al.|2005 ); a series of randomly 
orientated accretion events leads to a low spin value: we ex- 
pect an average value |a*| ~ 0.1-0.3 (King & Pringle 2006). 



The MBH also grows through mergers ( |Yu Sz Tremaine 2002 
Malbon et al.|2007| ). Minor mergers with smaller black holes 



(BHs) can decrease the spin (iHughes & Blandford||2003|) 
while a series of major mergers, between similar mass MBHs 
would lead to a likely spin of |a*| ~ 0.69 (Berti et al.||2007 



Gonzalez et al. 2007). Measuring the spin of MBHs shall 
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help us understand the relative importance of these pro- 
cesses, and perhaps gain a glimpse into their host galaxies' 
pasts. 

Elliptical and spiral galaxies are believed to host MBHs 
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of differing spins because of their different evolutions: we ex- 
pect MBHs in elliptical galaxies to have on average higher 
spins than black holes in spiral galaxies, where random, 
small accretion episodes have played a more important role 
(Volonteri, Sikora, &; Lasota 2007; Sikora, Stawarz, & Lasota 



2007 ). 

It has been suggested that the spin of the Galaxy's MBH 
could be inferred from careful observation of the orbits of 



stars within a few milliparsecs of the GC (Merritt et al 



|2010| , although this is complicated because of perturbations 
due to other stars, or from observations of quasi-periodic 
oscillations in the luminosity of flares believed to originate 
from material orbiting close to the innermost stable orbits 
jGenzel et al.||2003| |Hamaus et al.||2009|), though th ere are 
difficulties in interpreting these results ( |Psa ltis 2008). 

This latter method, combined with a disc-seismology 
model, has produced a value of the dimensionless spin of 
a* = 0.44 ± 0.08. To obtain this result |Kato et al.] ( |2010 ) 
have combined their observations of Sgr A* with observa- 
tions of galactic X-ray sources containing solar mass BHs, 
to find a best-fit unique spin parameter for all BHs. It is 
not clear that all BHs should share the same spin parame- 
ter; especially as the BHs considered here differ in mass by 
six orders of magnitude. Even if BH spin is determined by 
a universal process, we expect some distribution of spin pa- 
rameters ( King, Pringle, &; Hofmann 2008 ; Berti & Volonteri 



2008). Thus we cannot precisely determine the spin of the 
galactic centre's MBH from an average. 

The spins of MBHs in active galactic nuclei have been 
inferred using X-ray observations of Fe K emission lines 



(Miller 2007 McClintock et al.|2011 


) . This has been done for 


a handful of other galaxies' MBHs ( 


Brenneman & Reynolds 


2006;|Miniutti et al.||2009| |Schmoll et al.||2009| |de la Calle 


Perez et al.| 2010| Zoghbi et al.||2010| |Nardini et al.| 2011| 



to the maximal value for an extremal Kerr black hole. Typ- 
ical results are in the intermediate range of a* ~ 0.7 with 
an uncertainty of about 10% on each measurement. 

While we can use the spin of other BHs as a prior, to 
inform us of what we should expect for the Galaxy's MBH, 
it is desirable to have an independent observation, a direct 
measurement. 

An exciting means of inferring information about the 
MBH is through gravitational waves (GWs) emitted when 
compact objects (COs), such as stellar mass BHs, neutron 
stars (NSs), white dwarfs (WDs) or low mass main sequence 
(MS) stars, pass close by ( Sathyaprakash & Schutz 2009). 
A space-borne detector, such as the Laser Interferometer 
Space Antenna (LISA) or the evolved Laser Interferometer 
Space Antenna (eLISA), is designed to be able to detect 
GWs in the frequency range of interest for these encounters 
(Bender et al. 1998; Danzmann & Riidiger 2003; Jennrich 



dige 

et al. 20111 Amaro-Seoane et al.||2012t | T 1 Tn e identification 
of waves requires a set of accurate waveform templates cov- 
ering parameter space. Much work has already been done 
on the waveforms generated when companion objects inspi- 
ral towards an MBH ( |Glampedakisl|2005| |Barack||2009| ); as 



they orbit, GWs carry away energy and angular momen- 
tum, causing the orbit to shrink until eventually the CO 
plunges into the MBH. The initial orbits may be highly el- 
liptical and a burst of radiation is emitted during each close 
encounter. These are extreme mass-ratio bursts (EMRBs; 
Rubbo, Holley-Bockelmann, & Finn 2006). Assuming the 



companion is not scattered, and does not plunge straight 
into the MBH, its orbit evolves, becoming more circular, 
and it begins to continuously emit significant gravitational 
radiation in the LISA/ eLISA frequency range. The result- 



ing signals are extreme mass-ratio inspirals (EMRIs; Amaro- 
|Seoane et al.||2007| |Amaro-Seoane'||2012| ). 



Studies of these systems have usually focused upon the 
phase when the orbit is close to plunge and completes a large 
number of cycles in the detector's frequency band, allowing 
a high signal-to-noise ratio (SNR) to be accumulated. Here, 
we investigate high eccentricity orbits. These are the initial 
bursting orbits from which an EMRI may evolve, and are the 
consequence of scattering from two body encounters. The 
event rate for the detection of such EMRBs with LISA ha s 
been estimated to be as high as 15 yr _1 (Rubbo et al. 2006 ), 



although this has been sub sequently revised downwa rds to 
the order of 1 yr _1 (|Hopman, Freitag, fc Larson 2007). Even 



if only a single burst is detected during a mission, this is still 
an exciting possibility since the information carried by the 
GW gives an unparalleled probe of the spacetime of the GC. 
Exactly what can be inferred depends upon the orbit, which 
we investigate here. 

We make the simplifying assumption that all these or- 
bits are marginally bound, or parabolic, since highly eccen- 
tric orbits appear almost indistinguishable from an appropri- 
ate parabolic orbit. Here "parabolic" and "eccentricity" refer 
to the energy of the geodesic and not to the geometric shape 
of the orbit Following such a trajectory an object may 
make just one pass of the MBH or, if the periapsis distance 
is small enough, it may complete a number of rotations. 



Such an orbit is referred to as zoom-whirl (Glampedakis & 
Kennenck|2002 ). 



In order to compute the gravitational waveform pro- 
duced in such a case, we integrate the geodesic equations 
for a parabolic orbit in Kerr spacetime. We assume the or- 
biting body is a test particle, such that it does not influence 
the underlying spacetime, and that the orbital parameters 
evolve negligibly during the orbit such that they may be 
held constant. We use this to construct an approximate nu- 
merical kludge (NK) waveform ( Babak et al.||2007| ). 

This paper is organised as follows. We begin in Sec. [2] 
with the construction of the geodesic orbits; these trajecto- 
ries are used for the NK waveforms as explained in Sec. [3] 
In Sec. E] we establish what the LISA detectors would mea- 
sure and how the signal would be analysed. This includes 
a brief mention of window functions which is expanded in 
Appendix [A] Here we present a novel window function, the 
Planck-Bessel window, of use for signals with a large dy- 
namic range. In Sec. [5] we look at our NK waveforms. We 
give fiducial power-law fits for SNR as a function of peri- 
apse radius, useful for back-of-the-envelope estimates. We 
confirm the accuracy of the kludge waveforms in Sec. [6] by 



1 The revised eLISA concept shares the same descoped design as 
the New Gravitational- wave Observatory (NGO) submitted to 
the European Space Agency for their LI mission selection. 



2 Marginally bound Keplerian orbits in flat spacetime are 
parabolic in both senses. 
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comparing the energy flux to fluxes calculated using other 
approaches. The typical error introduced by the NK approx- 
imation may be a few percent, but this worsens as the pe- 
riapsis approaches the last non-plunging orbit. We explain 
how to extract the information from the bursts in Sec. 13 
Results estimating the measurement precision are presented 
in Sec. [8] We briefly mention the possibility of detecting 
bursts from extra-galactic sources in Sec. [5] before conclud- 
ing in Sec. [lO] with a summary of our results. EMRBs may 
be informative if the event rate is high enough for them to 
be a viable source. 

There are currently no funded space-borne detector mis- 
sions. The eLISA mission concept remains an active field of 
study. It is hoped to submit this to the European Space 
Agency as a potential cornerstone mission. We use the clas- 
sic LISA design. This is done from historical affection in lieu 
of a definite alternative. Should funding for a space-borne 
detector be secured in the future it is hoped that it shall 
have comparable sensitivity to LISA, and that studies using 
the LISA design shall be a sensible benchmark for compari- 
son. We find that to obtain good results the periapse radius 
must be r p < 10r g , where r g = GM m /c 2 is a gravitational 
radius; at this point the SNR is already high: for parame- 
ter estimation the orbit is more important that the signal 
strength, and so the exact detector performance should be 
of secondary importance. 

We adopt a metric with signature (+,—,—,—). Greek 
indices are used to represent spacetime indices \i — 
{0, 1, 2, 3} and lowercase Latin indices from the middle of the 
alphabet are used for spatial indices i = {1, 2, 3}. Uppercase 
Latin indices from the beginning of the alphabet are used 
for the output of the two LISA detector-arms A = {I, II}, 
and lowercase Latin indices from the beginning of the alpha- 
bet are used for parameter space. Summation over repeated 
indices is assumed unless explicitly noted otherwise. Geo- 
metric units with G = c = 1 are used where noted, but in 
general factors of G and c are retained. 



PARABOLIC ORBITS IN KERR 
SPACETIME 



here Loo is the total specific angular momentum at infinity, 



Astrophysical BHs are described by the Kerr metric (Kerr 
1963). This is conveniently written in terms of Boyer- 



Lindquist coordinates {t, r, fl, (j)} (|Boyer &; Lindquist 1967 
Hobson, Efstathio u, &; Lasenby|2006[ section 13.7). For this 



section we shall work in natural units with G = c = 1. 

Geodesies are parametrized by three conserved quan- 
tities (aside from the particle's mass /x): energy (per unit 
mass) E, specific angular momentum about the symme- 
try axis (the z-axis) L z , and Carter constant Q ( |Carter| 
1968; Chandrasekhar 1998, section 62). For a parabolic or- 
bit E — 1; the particle is at rest at infinity. This simplifies 
the geodesic equations. It also allows us to give a simple 
interpretation for the Carter constant: this is defined as 



1980 



where the metric is asymptotically flat (de Felice 
This is as in Schwarzschild spacetime. 

2.1 Integration variables and turning points 

In integrating the geodesic equations, difficulties can arise 
because of the presence of turning points, when the sign of 
the r or geodesic equation changes. The radial turning 
points are at the periapsis r p and at infinity. We may locate 
the periapsis by finding the roots of 



V r 2M.r 
0. 



(L 2 Z + Q) r 2 + 2M. [(L z 



- a) 2 + i 



(5) 

This has three roots, which we denote {ri,r2,r p }; the peri- 
apsis r p is the largest real root. We do not find the apoapsis 
as a (fourth) root as we have removed it by taking E — 1 
before solving]^] 

We avoid the difficulties associated with the turning 
point by introducing angular variables that always increase 
with proper time (Drasco & Hughes 2004): inspired by Ke- 
plerian orbits 

r T - 5 (6) 

1 + e cos ip ' 

where e = 1 is the eccentricity and p = 2r p is the semilatus 
rectum. As ip covers its range from — ty to 7r, r traces out a 
complete orbit. The geodesic equation for ip is 

dip 



(r 2 + a 2 cos 2 ( 



+ 



rir 2 
2r D 



2r p - (n + r 2 ) (1 + cos^) 

1/2 



(i + cos ipy 



(7) 



for proper time r. Parametrizing an orbit by its periapsis 
and eccentricity has the additional benefit of allowing easier 



comparison with its flat-space equivalent ( |Gair, Kennefick7 
fe Larson|2005[ ). 

The motion is usually bounded, with 0o ^ ^ tt — 0o; 
in the event that L z = the particle follows a polar orbit 



and covers its full range (Wilkins 1972). Turning points 
are given by 



V e = Q-cot 2 0L 2 z = 0. 

Changing variable to £ 
cos 2 0o given by 

Q = Q_ 

L 2 ^ 



(8) 



cos 2 6, we have a maximum £o 



(9) 



See Fig. [l] for a geometrical visualization. Introducing a sec- 
ond angular variable (Drasco & Hughes 2004) 

£ = £ cos 2 x; (10) 

over one 2tt period of x, oscillates from its minimum to its 
maximum and back. The geodesic equation for x is 

dr 



(r 2 + a 2 cos 2 ( 



(ii) 



Q = L 2 q + cos 2 I 



a 2 (1 



E 2 ) + ^ 



where Lq is the (non-conserved) specific angular momentum 
in the ^-direction. For E — 1 



Q = L 2 e + cot 2 L 2 



( 3 ) 3 See |Rosquist, Bylund, fe Samuelsson] ( |2009| ) for a discussion of 
the interpretation of Q in the limit G —> 0, corresponding to a 
flat spacetime. 

4 This turning point can be found by setting the unconstrained 
expression for V r equal to zero, and then solving for E(r); taking 

(4) the limit r oo gives E 1 ( Wilkins|1972 >. 
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Figure 1. The angular momenta Loo, L z and \fQ define a right- 
angled triangle. The acute angles are #o, the extremal value of the 
polar angle, and l, the orbital inclination |Glampedaki s, Hughes, | 
|fc Kennenck|2002| >. 



3 WAVEFORM CONSTRUCTION 

We can now calculate the geodesic trajectory. The orbit- 
ing body is assumed to follow this track exactly; we ignore 
evolution due to the radiation of energy and angular mo- 
mentum, which should be negligible for EMRBs. From this 
trajectory we calculate the waveform using a semirelativis- 
tic approximation flRuffini & Sasaki||1981|): we assume the 
particle moves along the Kerr geodesic, but radiates as if 
it were in flat spacetime. This quick-and-dirty technique is 
known as a numerical kludge (NK), and has been shown to 
approximate well results computed by more accurate meth- 
ods (Babak et al. 2007). It is often compared to a bead 



travelling along a wire. The shape of the wire is set by the 
geodesic, but the bead moves along in flat space. 



3.1 Kludge approximation 

Numerical kludge approximations aim to encapsulate the 
main characteristics of a waveform by using the exact par- 
ticle trajectory (ignoring inaccuracies from radiative effects 
and from the particle's self- force), whilst saving on com- 
putational time by using approximate waveform generation 
techniques. 

We build an equivalent flat-space trajectory by identify- 
ing the Boyer-Lindquist coordinates with a set of flat-space 
coordinates. We consider two choices: 

(i) Identify the Boyer-Lindquist coordinates with flat- 
space spherical polars {rsL, #bl, bl} —> {r sp h, fl sp h, sp h}, 
then define Cartesian coordinates ( |Gair et al . 2005; Babak 
et al. 112007) 



(r sp h sin sp h cos (p sp h \ 
r sph sin 6> sp h sin sph . (12) 
r sp h cos # sp h / 

(ii) Identify the Boyer-Lindquist coordinates with flat- 
space oblate-spheroidal coordinates {tbl 3 $bl 3 0bl } — >• 
{^ob, #ob, </>ob} so that the Cartesian coordinates are 

(Vr b 2 + a 2 sin oh cos o b\ 
\/r oh 2 + a 2 sin oh sin (j) oh . (13) 
r G b cos #ob / 

These are appealing because in the limit that G —> 0, where 
the gravitating mass goes to zero, the Kerr metric in Boyer- 
Lindquist coordinates reduces to the Minkowski metric in 
oblate-spheroidal coordinates. 

The two coincide for a — »• or r — ► oo. 

There is no well motivated argument that either coordi- 
nate system must yield an accurate GW; their use is justified 



post facto by comparison with results obtained from more ac- 
curate, and computationally intensive, methods ( |Gair et aTj 
2005 Babak et al. 2007). The ambiguity in assigning flat- 



space coordinates reflects the inconsistency of the semirel- 
ativistic approximation: the geodesic trajectory was calcu- 
lated for the Kerr geometry; by moving to flat spacetime we 
lose the reason for its existence. This should not be regarded 
as a major problem; it is an artifact of the basic assumption 
that the shape of the trajectory is important for determin- 
ing the character of the radiation, but the curvature of the 
spacetime in the vicinity of the source is not. By binding the 
particle to the exact geodesic, we ensure that the waveform 
has spectral components at the correct frequencies, but by 
assuming flat spacetime for generation of GWs they shall 
not have the correct amplitudes. 

3.2 Quadrupole-octupole formula 

Now we have a flat-space particle trajectory Xp(r), we may 
apply a flat-space wave generation formula. We use the 
quadrupole-octupole formula to calculate the gravitational 



strain ( |Bekenstein1|1973| |Press||1977| |Yunes et alT]|2QQ8| ) 

2G 



h jk {t,x) = 



2mS ijk + mM' 



1 

/ t' = t-r/c 



(14) 



where an over-dot represents differentiation with respect to 
time t, t' is the retarded time, r = \x — xp\ is the ra- 
dial distance, n is the radial unit vector, and I jk , S tjk 
and M % ^ k are the mass quadrupole, current quadrupole and 
mass octupole respectively. This is correct for a slowly mov- 
ing source. It is the familiar quadrupo le formula (|Misner, 
|Thorne, h Wheeler|1973l section 36.10; [Hobson et al.|2006[ 
section 17.9), derived from linearized theory, plus the next 
order terms. 



4 SIGNAL DETECTION AND ANALYSIS 
4.1 The LISA detector 

The classic LISA design is a three arm, space-borne laser 



interferometer (Bender et al. 1998 



2003). The arms form an equilateral triangle that rotates as 
the system's centre of mass follows a circular, heliocentric 
orbit, trailing 20° behind the Earth. eLISA has a similar 
design, but has only two arms, which are shorter in length, 
and trails 9° behind the Earth ( |Jennrich et al.||2011 ). 

To describe the detector configuration, and to trans- 
form from the MBH coordinate system to those of the de- 
tector, we use three coordinate systems: those of the BH at 
the GC x\\ ecliptic coordinates centred at the solar system 
(SS) barycentre x@, and coordinates that co-rotate with the 
detector x d . The MBH's coordinate system and the SS coor- 
dinate system are depicted in Fig. [2] The mission geometry 
for LIS A I eLISA is shown in Fig. [3] We define the detec- 
tor coordinates such that the detector-arms lie in the xa-yd 
plane as in ICutlerl (119981). We have computed the waveforms 



in the MBH's coordinates, but it is simplest to describe the 
measured signal using the detector's coordinates. 

The strains measured in the three arms can be combined 
such that LISA behaves as a pair of 90° inte rferometers a t 
45° to each other, with signals scaled by V3/2 ( |Cutler 1998 ). 
We denote the two detectors as I and II and use vector 
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introducing Fourier transforms 



Figure 2. The relationship between the MBH's coordinate sys- 
tem x\ and the SS coordinate system Xq. The MBH's spin axis 
is aligned with the z»-axis. The orientation of the MBH's x- and 
y-axes is arbitrary. We choose x m to be orthogonal to the direction 
to the SS. 




Figure 3. The relationship between the det ector coordinates x l r 
and the ecliptic coordinates of the SS Xq (Bender et al. 1998 
|Jennrich et al.|2011| ). The detector inclination is a = 6U U . 



notation h(t) — (hi(t) , hn(t)) 
from both detectors. 



to represent signals 



4.2 Frequency domain formalism 

Having constructed the GW h(t) that shall be incident upon 
the detector, we may consider how to analyse the waveform 
and extract the information it contains. We briefly recap 
GW signal analysis, with application to LISA. A more com- 



plete discussion can be found in Finn (1992) and Cutler & 
Flanagan (19941). Adaption for eLISA requires a substitu- 
tion of the noise distribution, and the removal of the sum 
over data channels, since it would only have one. 

The measured strain s(t) is the combination of the sig- 
nal and the detector noise 



s{t) = h(t)+n(t)- 



(15) 



we assume the noise 71a(£) is stationary and Gaussian, and 
that noise in the two detectors is uncorrelated, but shares 
the same characterisation ( Cutler||1998| . 

The properties of the noise allow us to define a natural 
inner product and associated distance on the space of signals 



(Cutler fe Flanagan|| 1994|) 



(g\k) 



Jo 



Sn(f) 



df, 



(16) 



3(f) = *{g(t)} 



-f 

J — ( 



g(t)e 



2-Kift 



dt, 



(17) 



and S n (f) is the noise spectral density. The signal-to- noise 
ratio is approximately 

p[h] = (h\h) 1/2 . (18) 

The probability of a particular realization of noise n(t) — 
n (t) is 



p(n(t) — no(t)) oc exp 



(n |n ) 



(19) 



Thus, if the incident waveform is h(t), the probability of 
measuring signal s(t) is 

p(s(t)\h(t)) oc exp -^-(s-h\s-h) . (20) 



4.3 Noise curve 

LISA's noise has two sources: instrumental noise and confu- 
sion noise, primarily from white dwarf binaries. The latter 
may be divided into contributions from galactic and extra- 
galactic binaries. We use the noise model of Barack & Cutler 
(2004J. The instrumental noise dominates at both high and 
low frequencies. The confusion noise is important at inter- 
mediate frequencies, and is responsible for the cusp around 
10 -3 Hz. eLISA shares the same sources of noise, but is 
less affected by confusion. Its sensitivity regime is shifted to 
higher frequencies because of a shorter arm length. 

4.4 Window functions 

There is one remaining complication regarding signal anal- 
ysis: as we are Fourier transforming a finite signal we en- 
counter spectral leakage; a contribution from large ampli- 
tude spectral components leaks into surrounding compo- 
nents (sidelobes), obscuring and distorting the spectrum at 
these frequencies (Harris 1978). This is an inherent problem 
with finite signals; it shall be as much of a problem when 
analysing signals from an actual mission as it is here. To 
mitigate, but unfortunately not eliminate, these effects, the 
time-domain signal can be multiplied by a window function. 
These are discussed in detail in Appendix^] We adopt the 
Nuttall four-term window with continuous first derivative 



flNuttall|1981t for our results. 



5 WAVEFORMS AND DETECTABILITY 
5.1 Model parameters 

The shape of the waveform depends on parameters defining 
the MBH; the companion object on its orbit, and the de- 
tector. Let us define A = {A 1 , A 2 , . . . , A z } as the set of Z 
parameters which specify the GW. For our model, these are: 

(1) The MBH's mass M m . This is well constrained by 
the observation of stellar orbits about Sgr A* ( |Ghe: 
[2009t 



et al. 



^008^ IGillessen 



al. 



with best estimate M m — 
1X31 ±0.36) x 1O°M . Thisdepends upon the galactic cen- 
tre distance R as M. = (3.95 ± 0.06| s tat ± 0.18| Rq , stat 
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0.31|h , sys ) x 10 6 M o (i? /8 kpc) 219 , where the errors are 
statistical, independent of Ro; statistical from the determi- 
nation of i?o, and systematic from Ro respectively. 

(2) The spin parameter a*. Naively this could be any- 
where in the range |a*| < 1; however, it is possible to 
place an upper bound by contemplating spin-up mecha- 
nisms. Considering the torque from radiation emitted by an 
accretion disc, and swallowed by the BH, it can be shown 
that |a*| < 0.998 flThorne|[l974|) . Magnetohydrodynamical 



simulations of accretion discs produce a smaller maximum 
value of | a* | - 0.95 ( jGammie, Shapiro, & McKinney|2004| ). 
The actual spin value could be much lower than this upper 
bound,? depending upon the MBH's evolution. 

(3, 4) The orientation angles for the MBH 9k and <3>k- 

(5) The ratio of the SS-GC distance Ro and the CO mass 
/i, which we denote as £ = Ro/ P- This scales the amplitude 
of the waveform. Bursts, unlike inspirals, do not undergo 
orbital evolution, hence we cannot break the degeneracy in 
Ro and p. The distance, like M # , is constrained by stellar 
orbits, the best estimate being (Gillessen et al. 2009) Ro — 
8.33 =b 0.35 kpc. The mass of the orbiting particle depends 
upon the type of object: whether it is an MS star, WD, NS 
or BH. Since we shall not know the p precisely, we shall not 
be able to infer anything more about the distance. 

(6, 7) The angular momentum of the CO. This can be 
described using either {L Z ,Q} or {Loo,^}. We employ the 
latter, as the total angular momentum and inclination are 
less tightly correlated. Assuming spherical symmetry, we ex- 
pect cos l to be uniformly distributed. 

(8-10) A set of coordinates to specify the trajectory. We 
use the angular phases at periapse, (j) p and Xp (which deter- 
mines P ), as well as the time of periapse t p . 

(11, 12) The coordinates of the MBH from the SS 
bary centre O and These may be taken as the coordi- 
nates of Sgr A*, as the radio source is expected to be within 
20r g of the MBH ( |Reid et al.||2003| |Doeleman et al .12008). 



We use the J2000.0 coordinates (Reid et al. 1999 lYusef- 



Zadeh et al.||1999 ). These change with time due to the ro- 



tation of the SS about the GC; the proper motio n is about 
6 masyr -1 , mostly in the plane of the galaxy (Backer & 
Sramek 1999; Reid et al. 2003). The position is already de- 
termined to high accuracy: an EMRB can only give weak 
constraints on source position^] We take it as known and do 
not try to infer it. 

(13, 14) The orbital position of the detector given by 
4> and ip. We assume the i nitial position s are chosen such 
that = when cp = (Cutler 1998); this choice does 



not qualitatively influence our results. The orbital position 
should be known, so need not be inferred. 



We therefore have a 14 dimensional parameter space, of 
which we are interested in inferring d = 10 parameters. 



5.2 Waveforms 

Figure [4] shows example waveforms to demonstrate some 
possible variations in the signal. All assume the standard 



5 For comparison, an EMRI, which should be more informative, 
can only give sky localisation to ~ 10 _ 3 steradians jBarack "&] 
Cutler|2004| [Huerta fc Gair|2009| >. 



MBH mass and position as well as a p — 10 Mq CO; other 
(randomly chosen) orbital parameters are specified in the 
captions. Radii are given in gravitational radii r g = GM m /c 2 . 

The plotted waveforms use the spherical polar coordi- 
nate system for the NK. Using oblate-spheroidal coordinates 
makes a small difference: on the scale shown the only dis- 



cernible difference would be in Fig. 4(c) the maximum dif- 
ference in the waveform (outside the high-frequency tail) is 
~ 10%. In the other cases the difference is entirely negligible 
(except in the high-frequency tail, which is not of physical 
significance). This behaviour is typical; for the closest or- 
bits, with the most extreme spin parameters, the maximum 
difference in the waveforms may be ~ 30%. The difference is 
largely confined to the higher frequency components, which 
are most sensitive to the parts of the trajectory closer to the 
MBH: the change in flat-space radius for the same Boyer- 
Lindquist radial coordinate causes a slight shift in the shape 
of the spectrum. Enforcing the same flat-space periapse ra- 
dius gives worse agreement across the spectrum. 

To examine the effect of the coordinate choice, we com- 
pare SNRs calculated using the alternative schemes. The 
MBH parameters were fixed as for the GC, the orbital pa- 
rameters were chosen such that periapse distance was drawn 
from a logarithmic distribution (down to the innermost sta- 
ble orbit), and other parameters were drawn from appro- 
priate uniform distributions. The ratio of the two SNRs is 
shown in Fig. [5] The difference from the coordinate systems 
is only apparent for orbits with very small periapses. There 
is agreement to 10% down to r p ~ 4r g ; the maximal differ- 
ence may be expected to be ~ 20%, this is for periapses that 
are only obtainable for high spin values. 

Since the deviation in the two waveforms is only ap- 
parent for small periapses, when the kludge approximation 
is least applicable, we conclude that the choice of coordi- 
nates is unimportant. The potential error of order 10% is 
no greater than that inherent in the NK approximation (see 
Sec.[6|. Without an accurate waveform template to compare 
against, we do not know if there is a preferable choice of 
coordinates. We adopt spherical coordinates for easier com- 
parison with existing work. 

5.3 Signal-to-noise ratios 

The detectability of a burst depends upon its SNR. To char- 
acterise the variation of p we considered a range of orbits, 
each for the standard MBH mass and position. These bursts 
were calculated for a 1M© CO. From (14), the amplitude of 



the waveform is proportional to the CO mass p and so p is 
also proportional to p; a 10M© object is ten times louder on 
the same orbit. To make results mass independent, we use 
a mass-normalised SNR 



p[h} = 



M 



p[h]. 



(21) 



The spin of the MBH and the orbital inclination were 
randomly chosen, and the periapse distance was drawn from 
a logarithmic distribution down the inner-most stable orbit. 
For each set of these extrinsic parameters, the periapse posi- 
tion, orientation of the MBH, and orbital position of the de- 
tector were varied: five random combinations of these intrin- 
sic parameters (each being drawn from a separate uniform 
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(a) Waveform for a* ~ 0.12, r p c± 15.6r g and i ~ 2.1. The 
SNR for the spherical polar kludge is p[fa sp h] — 451, for 
the oblate-spheroidal kludge p[fa b] — 451 (agreement to 
0.01%). 
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Figure 5. Ratio of SNR for a waveform calculated using spherical 
polar coordinates to that for a waveform using oblate-spheroidal 
coordinates. 
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(b) Waveform for a* ~ 0.48, r p ~ 8.8r g and l ~ 2.0. The 
SNR for the spherical polar kludge is p[^i S ph] — 2300, for 
the oblate-spheroidal kludge p[fo b] — 2310. 
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(c) Waveform for a* ~ 0.74, r p ~ 3.2r g and i ~ 1.2. The 
SNR for the spherical polar waveform is p[fa sp h] — 70600, 
for the oblate-spheroidal kludge p[fo Q b] — 74900. 

Figure 4. Example burst waveforms from the galactic centre. 
The strain h\(f) is indicated by the solid line, hn(f) by the dot- 
dashed line, and the noise curve by the dashed line. The kludge 
has been formulated using spherical polar coordinates. 
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r p/ r g 



Figure 6. Mass-normalised SNR as a function of periapse ra- 
dius. The plotted points are the values obtained by averag- 
ing over each set of intrinsic parameters. The best fit line is 
log(p) = -2.691og(r p /r g ) +4.88. This is fitted to orbits with 
r p > 13.0r g and has a reduced chi-squared value of x 2 l u — 1-73. 

distribution) were used for each point. We take the mean of 
In p for each set of randomised intrinsic parameters^] 

There exists a correlation between the periapse radius 
and SNR, as shown in Fig. [6] Closer orbits produce louder 
bursts. To reflect this trend, we have fitted a simple fiducial 
power law, indicated by the straight lineQThis was done by 
maximising the likelihood, assuming In p has a Gaussian dis- 
tribution with standard deviation derived from the scatter 
because of variation in the intrinsic parameters. The power 
law is a good fit only for larger periapses. The shape is pre- 
dominately determined by the noise curve. The change in the 

6 The logarithm is a better quantity to work with since the SNR 
is a positive-definite quantity that may be distributed over a range 
of magnitudes (MacKay 2003 sections 22.1, 23.3). Using median 
values yields results that are quantitatively similar. 

7 Using oblate-spheroidal coordinates instead of spherical polars 
gives a fit consistent to within 0.1% as we have excluded the 
closest orbits. 
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trend reflects the transition as from approximately power 
law behaviour to the bucket of the noise curve. Hence, we 
fit a power law to orbits with a characteristic frequency of 
/* = ^/GM m /r p < 1 x 1CT 3 Hz, to avoid spilling into the 
bucket. Changing the cut-off within a plausible region alters 
the fit coefficients by around 0.l[^] 

The SNR shows no clear correlation with the other pa- 
rameters (excluding p). However, the SNR is sensitive to the 
intrinsic parameters, in particular the initial position, and 
may vary by an order of magnitude. 

Setting a threshold of p = 10, a 1M (10M©) object 
would be expected to be detectable if the periapse distance 
is less than 27r g (65r g ). Hopman et al. (2007), assuming a 



threshold of p = 5, used an approximate form for the SNR 
based upon the quadrupole component of a circular orbit; 
their model, with updated parameters for the MBH, predicts 
bursts would be detectable out to 66r g (135r g ). This is overly 
optimistic. 



6 ENERGY SPECTRA 

To check the NK waveforms, we compare the energy spectra 
calculated from these with those obtained from the classic 

(fT964|. 



treatment of Peters & Mathews (1963) and 



Peters 



This calculates GW emission for Keplerian orbits in flat 
spacetime, assuming only quadrupole radiation. The spec- 
trum produced should be similar to that obtained from the 
NK in weak fields, that is for large periapses; we do not ex- 
pect an exact match because of the differing input physics 
and varying approximations. 

In addition to using the energy spectrum, we also 
use the total energy flux. This contains less information 
than the spectrum; however, results have been calculated 
for parabolic orbits in Schwarzschild spacetime using time- 
domain black hole perturbation theory (Martel 2004). These 
should be more accurate than results calculated using the 
Peters and Mathews formalism. 

We do not intend to use the kludge waveforms to cal- 
culate an accurate energy flux: this would be inconsistent as 
we assume the orbits do not evolve with time. We only cal- 
culate the energy flux as a sanity check, to confirm that the 
kludge approximation is consistent with other approaches. 



6.1 Kludge spectrum 

A gravitational wave in the TT gauge has an effective 
energy-momentum tensor ( Misner et al.|1973| section 35.15) 



T 



32ttG 



(22) 



where (...) indicates averaging over several wavelengths or 
periods. The energy flux through a sphere of radius R is 



dE 
dt 



32nG 



R 



I 



dt dt 



(23) 



8 The power law exponent —2.7 is inconsistent with —13/4 as 
predicted by the approximate model of Hopma n et al.| ( |2007| ). 
This is the result of their approximate waveform model. 



with f dQ representing integration over all solid angles. 
From |l4| the waves have a 1/r dependence; if we define 



hi 



Hi 



(24) 



we see, using (14), the flux is independent of i?, as required 
for energy conservation, and 



dE _ c 
dt ~ 32 



nGj 



dn/ d^d^ 



dt dt 

Integrating to find the total energy emitted we obtain 



E 



32nG 



dt , : . 

dt dt 



(25) 



(26) 



Since we are considering all time, the localization of the 
energy is no longer of importance and it is unnecessary to 
averagejDver several periods. Switching to Fourier represen- 
tation Hij(f) = &{Hij(t)}, 



4G 



I dn L° 



d// 2 ir j ' (/)#*(/), 



(27) 



using H*j(f) — Hij(-f) as the signal is real. From this we 
identify the energy spectrum as 

3 J dn/ 2 # y (/)#* (/). (28) 



dE _ 7TC- 

df ~ 4G 



6.2 Peters and Mathews spectrum 

To calculate the Peters and Mathews energy spectrum for a 



parabolic orbit, we use the limiting result of Turner (1977). 
This result should be accurate to ~ 10% for orbits with 
periapse radii larger than ~ 20r g ( Berry &; Gair|20T0 ). 



6.3 Comparison 

Two energy spectra are plotted in Fig. [7] for orbits with 
periapses of r p = 15.0r g , 30.0r g and 60.0r g . The two spectra 
appear to be in good agreement, showing the same general 
shape in the weak- field limit. The NK spectrum is more 
tightly peaked, but is always within a factor of 2 at the 
apex. The peak of the spectrum is shifted to a marginally 
higher frequency in the NK spectrum primarily because of 
the addition of the higher order terms. 

Comparing the total energy fluxes, ratios of the various 
energies are plotted in Fig. [8] We introduce an additional 
energy, the quadrupole NK energy ^nk(Q)- This allows eas- 
ier comparison with the Peters and Mathews energy which 
includes only quadrupole radiation. It can be calculated in 
three ways: 

(i) Inserting the waveform h(f) generated including only 
the mass quadrupole term in |l4| into ( 27 ) and integrating. 



This is equivalent to the method used to calculate Enk- 

(ii) Numerically integrating the quadrupole GW luminos- 
ity ( |Misner et al. 1973| section 36.7; Hobson et al. 2006 
section 18.7) 



G 



I 



where Ii 



dt, 



(29) 



(l/3)ISij is the reduced mass quadrupole 



tensor. We can obtain this from (26), by integrating over 
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Figure 7. Energy spectra for a parabolic orbit of a 11 = IOMq object about a Schwarzschild BH with M m = 4.31 x 1O 6 M0. The spectra 
calculated from the NK waveform are shown by the solid line and the Peters and Mathews flux is indicated by the dashed line. The NK 
waveform includes current quadrupole and mass octupole contributions. The high frequency tail is the result of spectral leakage. 



all angles when the waveform only contains the mass 
quadrupole component. This has the advantage of avoid- 
ing the effects of spectral leakage or the influence of window 
functions. 



(iii) Using the analytic expressions for the integral (29) 
given in appendix A of Gair et al. (2005). 



All three agree to within computational error. No difference 
is visible on the scale plotted in Fig. [8] This demonstrates 
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r p /r g 2tt/A0 
(a) Versus periapsis (b) Versus amount of rotation 

Figure 8. Ratios of energies as a function of periapsis r p and 2tt divided by the total angle of rotation in one orbit Acf) (27r/A</> = 1 
for a Keplerian orbit). The solid line shows the ratio of the numerical kludge and Martel energies the dashed line shows the 

ratio of the NK energy calculated using only the mass quadrupole term and the Martel energy £^ NK (q)/^m; the dot-dashed line shows 
the ratio of the quadrupole and quadrupole-octupole NK energies ^nk(Q)/-^nk, and the dotted line shows the ratio of the Peters and 
Mathews and quadrupole NK energies ^pm/-Snk(Q)- The spots show the mapping between the two abscissa scales. Compare with figure 
4 of |Gair et al"| ( [2005| >. 



the validity of the code, and shows that the use of a window 
function does not significantly distort the waveform. 

The ratios all tend towards one in the weak field, as 
required, but differences become more pronounced in the 
strong field. The NK energy is larger than the Peters and 
Mathews result Epm- This behaviour has been seen before 



for high eccentricity orbits about a non-spinning BH (Gair 
et al. 2005). It may be explained by considering the total 



path length for the different orbits: the Peters and Mathews 
spectrum assumes a Keplerian orbit, the orbit in Kerr geom- 
etry rotates more than this. The greater path length leads 
to increased emission of gravitational waves and a larger en- 
ergy flux ( Ber ry &; Gair|2010 ) . Our bead must travel further 
along its wire. A good proxy for the path length is the an- 
gle of rotation A0; this is 2tt for a Keplerian orbit, in Kerr 
the angle should be 2tt in the limit of an infinite periapsis, 
whereas for a periapsis small enough that the orbit shows 
zoom-whirl behaviour, the total angle may be many times 
2tt. There is a reasonable correlation between the amount of 
rotation 2tt/A0 and the ratio of energies. 

Error in the NK energy compared with the time-domain 
black hole perturbation theory results of Martel comes from 
two sources: the neglecting of higher order multipole con- 
tributions and the ignoring of background curvature. The 
contribution of the former can be estimated by looking at 
the difference in the NK energy by including the current 
quadrupole and mass octupole terms. From Fig. [8] we see 
these terms give a negligible contribution in the weak field, 
but the difference is ~ 20% in the strong field. This ex- 
plains why the Martel energy Em is greater in the strong 
field: it includes contributions from all multipoles. Neglect- 
ing the background curvature increases the NK energy rel- 
ative to Em- This partially cancels out the error introduced 
by not including higher order terms: this accidentally leads 
to E NK (q) being more accurate than Enk for r p > 10r g 
( Tanaka et al.|1993| ). 



From the level of agreement we may be confident that 
the NK waveforms are a reasonable approximation. The dif- 
ference in energy flux is only greater than 10% for very 
strong fields r p ~ 4r g ; since this is dependent on the square 
of the waveform, typical accuracy in the waveform may be 
- 5% ( |Gair et aI.|2005||Tanaka et al.|1993| ). This is more sig- 
nificant than the variation in waveforms we generally found 
using the two alternative coordinate systems for the NK (in 
this case the two coincide because a* =0). 



7 PARAMETER ESTIMATION 

Having detected a signal, we are interested in what we can 
learn about the source. We have an inference problem that 



can be solved by application of Bayes' Theorem (Jaynes 



2003 chapter 4): the probability distribution for our param- 



eters given that we have detected the signal s(t) is given by 
the posterior 

p(s(t)\\)p(\) 



P (\\s(t)) 



p(s(t)) 



(30) 



Here p(s(t)\X) is the likelihood of the parameters, p(\) is 
the prior probability distribution for the parameters, and 
the evidence p(s(t)) — f p(s(t)\X) d d X is, for our purposes, 
a normalising constant. The likelihood depends upon the 
realization of noise. If parameters Ao define a waveform 
ho(t) = h(t; Ao), the probability that we observe signal s(t) 



GW is given by ( 20 ) , so the likelihood is 



p(s(t)\X ) oc exp 



(s - h \s - h ) 



(31) 



If we were to define this as a probability distribution for the 
parameters A, the modal values are the maximum-likelihood 
(ML) parameters Aml- The waveform h(t; Aml) is the signal 
closest to s(£), where distance is denned using the inner 
product (16) ( Cutler fe Flanagan||1994 ). 
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7.1 Fisher matrices 



and Ah is the corresponding change in the waveform 



In the limit of a high SNR, we may approximate (Vallisneri 
20081 



Ah h(\i a ) - h(Xb 



(39) 



p(s(t)\\ ) oc exp - 1 - (d a h\d b h) (X a - (\ a )z) (A 6 - (A b ) J 

(32) 

where the mean is defined as 
f X a p(s(t)\X) d d X 



fp{a(t)\\)d*\ ' 



(33) 



In the high SNR limit, this is the ML value (X a ) e = Aml- 
The quantity 



r ab = {d a h\d b h) 



(34) 



is the Fisher information matrix (FIM). It controls the vari- 
ance of the likelihood distribution. 

The form of the posterior distribution depends upon the 
nature of the prior information. If we have an uninformative 
prior, such that p(A) is a constant, the posterior distribution 
is determined by the likelihood. In the high SNR limit, we 
obtain a Gaussian with variance-covariance matrix 



s = r 1 . 



(35) 



The FIM therefore gives the uncertainty associated with the 
inferred parameters, in this case the ML values. 

If the prior restricts the allowed range for a parameter, 
as is the case for the spin a* , then the posterior is a truncated 
Gaussian, and T _1 may no longer represent the variance- 
covariance. 

If the prior is approximately Gaussian with variance- 
covariance matrix So, the posterior is also Gaussian ^] The 



posterior variance-covariance is (Cutler & Flanagan 
Vallisneri||2008| ) 



1994) 



e= (r + Eo- 1 ) 



(36) 



From this the inverse FIM T _1 is an upper bound on the 
size of the posterior covariance matrixp^] 

The FIM gives a quick way of estimating the range of 
the posterior. It is widely used because of this. However, it 
is only appropriate when the approximation of ( [32] ) holds. 
This is known as the linearised-signal approximation (LSA), 
where higher order derivatives are neglected. To assess the 
validity of this, |Vallisneri| ( |2008| recommends use of the 
maximum-mismatch (MM) criterion 

lnr = -i (A\ a d a h M L - Ah\A\ b d b h M L - Ah}) . (37) 
Here A A is the displacement to some point on the la surface 



AA = Ai 



(38) 



9 If we only know the typical value and spread of a parameter, 
a Gaussian is the maximum entropy prior (Jaynes 2003] section 
7.11): the prior that is least informative given what we know. 

10 It is also the Cramer-Rao bound on the error covariance of 
an unbiased estimator (Cutler & Flanagan 1994 Vallisneri 2008). 
Thus it represents the frequentist error: the lower bound on the 
covariance for an unbiased parameter estimator A es t calculated 
from an infinite set of experiments with the same signal h(t) but 
different realisations of the noise n(t). 



The la surface is defined from the inverse of the FIM. If 
higher order terms are indeed negligible, the MM criterion is 
small. We check this by picking a random selection of points 
on the la surface and evaluating | lnr|. If this is smaller than 
a fiducial value (| lnr| = 0.1) over the majority (90%) of the 
surface we consider the LSA sufficiently justified. 

We calculated FIMs for a wide range of orbits and 
checked the MM criterion. We found that for the overwhelm- 
ing majority the test failed: the LSA is not appropriate. 
This behaviour was seen even for orbits with p ~ 10 3 -10 4 p^| 
Higher order terms are important, and cannot be neglected. 

EMRBs have a short duration and accordingly are not 
the most informative of signals. Therefore, the la surface as 
defined by considering only the LSA terms is large. Taking 
such a step in parameter space moves the signal beyond the 
region of linear changes. 

We hope that this shall serve as an example to others. 
What constitutes high SNR depends upon the signal; it is 
not enough for p > 1. As stressed by Vallisneri (2008), it is 



essential to check the MM criterion for individual waveforms: 
the threshold for the LSA to become applicable could be 
much greater than naively thought. 

As we cannot be confident in FIM results, we aban- 
don this approach in favour of using Markov chain Monte 
Carlo simulations to explore constraints from different re- 
gions of parameter space. These are computationally more 
expensive, but do not rely on any approximations. 

7.2 Markov chain Monte Carlo methods 

Markov chain Monte Carlo (MCMC) methods are widely 
used for inference problems; they are a family of algorithms 
for integrating over complicated distributions and are effi- 
cient for high-dimensional problems (MacKay 2003, chapter 
29). Parameter space is explored by constructing a chain of 
N samples. The distribution of points visited by the chain 
maps out the underlying distribution; this becomes asymp- 
totically exact as N — >• oo. Samples are added sequentially, if 
the current state is A n a new point A* is drawn and accepted 
with probability 



A 



l 1 

{ 7TI 



(A*)£(A*)Q(A n ; A* 
(A n )£(A n )Q(A n ; A* 

setting A n +i = A*, where C(X) 
case from (31); 7r(A) is the prior, 



i4 



(40) 



is the likelihood, in our 
and Q is a proposal dis- 
tribution. If the move is not accepted A n +i = A n . This is 



the Metropolis-Hastings algorithm (Metropolis et al.||1953 



Hastings|1970) 



Waiting long enough yields an exact posterior, but it is 
desirable for the MCMC to converge quickly. This requires 
a suitable choice for the proposal distribution, which can be 
difficult, since we do not yet know the shape of the target 
distribution. 

11 In this study, to increase p we must reduce the periapse dis- 
tance; this also reduces the region where the LSA is valid as pa- 
rameter dependencies become more non-linear. If we had the lux- 
ury of increasing p by moving the GC closer, things could be 
different. 
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One method to define the proposal is to use the previ- 
ous results in the chain and refine Q by learning from these. 
Such approaches are known as adaptive methods. Updat- 
ing using previous points means that the chain is no longer 
Markovian. Care must be taken to ensure that ergodicity is 
preserved and convergence obtained (iRoberts & Rosenthal 



2007 Andrieu & Thorns 2008). To avoid this complication, 



we follow Haario et al. (1999), and use the adapting method 



as a burn in phase. We have an initial phase where the pro- 
posal is updated based upon accepted points. After this we 
fix the proposal and proceed as for a standard MCMC. By 
only using samples from the second part, we guarantee that 
the chain is Markovian and ergodic, whilst still enjoying the 
benefits of a tailor-made proposal. After only a finite number 



of samples we cannot assess the optimality of this ( [Andrieu 
&; Thoms|2008| ), but the method is still effective 



To tune Q, we use an approach based upon the adap- 



tive Metropolis (AM) algorithm ( |Haario et aL||2001| ). The 
proposal is taken to be a multivariate normal distribution 
centred upon the current point, the covariance of which is 



C = s(V n + £C ) ; 



(41) 



where V n is the covariance of the accepted points 
{Ai, . . . , A n }, s is a scaling factor that controls the step size, 
e is a small positive constant (typically 0.0025), and Co is 
a constant matrix included to ensure ergodicity. 

Our adaptation is run in three phases. The initial phase 
is to get the chain moving. For this C nit is a diagonal ma- 
trix with elements calibrated from initial one dimensional 
MCMCs. This finishes after iVmit accepted points. 

For the second phase, we use the proposal covariance 
from the initial phase C init for C™ ain . We reset the covari- 
ance of the accepted points so that it only includes points 
from this phase. This is the main adaptation phase and lasts 
until Amain points have been accepted. 

In the final adaptation phase we restart the chain at 
the true parameter values. We no longer update the shape 
of the covariance (V n remains fixed), but adjust the step 
size s to tune the acceptance rate; it is then fixed, along 
with everything else, for the final MCMC. 

Throughout the adaptation, we update the step size s 
after every 100 trial points (whether or not they are ac- 
cepted). While updating, the covariance V n changes af- 
ter every 1000 trial points. We set Ai n it = 50000 and 
N main = 450000. 

We initially aimed for an acceptance rate of 0.234; this is 
optimal for a random walk Metropolis algorithm with some 



specific high-dimensional target distributions ( |Roberts et al. 
1997| |Roberts fe Rosenthal||2001| ). In many cases we found 
better convergence when aiming for a lower acceptance rate, 
say 0.1. This is not unexpected: the optimal rate may be 
lower than 0.234 when the parameters are not independent 
and identically distributed (Bedard 2007, 2008b a). In prac- 
tice, the final acceptance rate is (almost always) lower than 
the target rate as the use of a multivariate Gaussian for the 
proposal distribution is rarely a good fit at the edges of the 
posterior. Consequently, the precise choice for the target ac- 
ceptance rate is unimportant as long as it is of the correct 
magnitude. Final rates are typically within a factor of 2 of 
the target value. As an initial choice, we set s = 2.38 2 /d, 
which is the optimal choice if C was the true target covari- 
ance for a high dimensional target of independent and iden- 



tically distributed parameters ([German et al.| 1996] IRob erts 
et al.|l997| [Roberts fe Rosenthal|2001||Haario et al.|200 1)p 
To assess the convergence of the MCMC we check the 
trace plot (the parameters values throughout the run) for 
proper mixing, that the one and two dimensional posterior 
plots fill out to a smooth distribution, and that the distri- 
bution widths tend towards consistent values. 



8 RESULTS 
8.1 Data set 

To investigate the information contained in EMRBs, we 
again considered a range of orbits. The MBH was assumed 
to have the standard mass and position. The CO was chosen 
to be 10M©, as the most promising candidates for EMRBs 
would be BHs: they are massive and hence produce higher 
SNR bursts, they are more likely to be on close orbits as 



a consequence of mass segregation (Bahcall &; Wolf 1977 



Alexander & Hopman| 2009i Pr eto k Amaro-Seoa ne 2010), 
and they cannot be tidally disrupted. 

Other potential source candidates would be NSs and 
low mass MS stars. Like BHs, NSs are not tidally disrupted 
in the GC, and benefit (albeit weakly) from mass segre- 
gation. MS stars have the great advantage of being more 
numerous than stellar remnants. Most MS stars would be 
tidally disrupted at periapses of interest, but the lowest 
mass stars c ould be viable sources (Freitag 2003; Amaro- 
Seoane| [2012 ) p] These could produce more intricate wave- 
forms, since they are extended bodies, subject to tidal dis- 
tortion, rather than point masses. We stick to black holes as 
the most credible and simplest source. Modelling extended 
bodies is left for future work. 

Orbits were chosen with periapses uniformly distributed 
in logarithmic space between the the inner-most orbit and 
16r g . The other parameters were chosen randomly from ap- 
propriate uniform distributions. 

The results of the MCMC runs illustrate why the FIM 
approach was insufficient. There are strong and complex pa- 
rameter dependencies. For many sets of parameters the pos- 
teriors are far from Gaussian as assumed in the LSA. Some 
example results are shown in Fig. [9] [To] and [XT] 

The first is well-behaved. It is almost Gaussian, but 
we see some asymmetries and imperfections. There are also 
strong degeneracies, indicated by needle-like distributions. 
This is a fairly standard example: there are runs which are 
closer to being Gaussian (especially at higher SNR), and 
equally there are tighter degeneracies. The lenticular M m - 
Loo degeneracy is common. 

The second shows banana-like degeneracies. These are 
not uncommon; there are varying degrees of curvature. The 
more complicated shape makes it harder for the MCMC to 
converge, so the final distribution is not as smooth as for 
the first example. The curving degeneracies also bias the one 
dimensional marginalisations away from the true values. 



12 Reasonably good results may be obtained by fixing s at this 
value, and not adjusting to fine tune the acceptance rate. 

13 The optimal density to resist disruption occurs at /i ~ O.O7M0 
(Chabrier & Baraffe 2000J. Bursts for these objects should be 
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Figure 9. Marginalised one and two dimensional posteriors. The scales are identical in both sets of plots. The dotted line indicates the 
true value. These distributions are fairly cromulent and well converged. Angular momentum is in units of L m = GM m c~ 1 . The input 
orbit has r p = 8.54r g and p = 916. 



The third shows more intricate behaviour. This is more 
rare, but indicates the variety of shapes that is obtainable. 
Again the convergence is more difficult, so the distributions 
are rougher around the edges; there is also some biasing due 
to the curving degeneracies. 

These results do not incorporate any priors (save to 
keep them within realistic ranges). We have not folded in 
the existing information we have, for example, about the 
MBH's mass. Therefore, the resulting distributions charac- 
terise what we could learn from EMRBs alone. By the time 
a space-borne GW detector finally flies, we will have much 
better constraints on some parameters. 

It is possible to place good constraints from the clos- 
est orbits. These can provide sufficient information to give 
beautifully behaved posteriors although significant correla- 
tion between parameters persists. 



8.2 Distribution widths 

Characteristic distribution widths are shown in Fig. |12[ 
Plotted are the standard deviation <tsd; a scaled 50- 
percent ile range a^o = W$q/1. 34898, where W50 is the range 
that contains the median 50% of points, and a scaled 95- 
percentile range (795 = W95/3. 919928, where W59 is the 
95% range. These widths are equal for a normal distribu- 
tion. Filled circles are used for runs that appear to have 
converged. Open circles are for those yet to converge, but 
which appear to be approaching an equilibrium state; widths 
should be accurate to within a factor of a few. For guidance, 
the dotted line corresponds to the current measurement un- 
certainty for M # ; the dashed lines are from uniform priors 
for a*, <J>k, P , Xp? cos Ok and cos^, and, for completeness, 
the solid line indicates the total prior range. We have no ex- 
pectations for the width of the MBH mass distribution with 
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Figure 10. Marginalised one and two dimensional posteriors. The scales are identical in both sets of plots. The dotted line indicates the 
true value. These distributions show definite non-gaussianity. The input orbit has r p = 9.86r g and p = 1790. 



respect to the current value; however, we would expect that 
the recovered distributions for the other parameters are nar- 
rower than for the case of complete ignorance. This may not 
be the case if the distribution is multimodal: in this event 
using the width is an inadequate description of the distribu- 
tion. Only a few unconverged runs exceed these limits, and 
some appear to be multimodal. 

The widths show a trend of decreasing with decreas- 
ing periapsis or increasing SNR, but there is a large degree 
of scatter. There does not appear to be a strong depen- 
dence upon any single input parameter, with the exception 
of the spin. The widths for l. Ok, ^k, </> p and %p increase 
for smaller spin magnitudes. The dependence is shown in 
Fig. [13] These parameters are defined with reference to the 
coordinate system established by the spin axis: for a* = 
we have spherical symmetry and there would be ambiguity 
in defining them. Therefore, it makes sense that they can 



be more accurately determined for larger spin magnitudes. 
The width for a*, however, shows no clear correlation. 

Comparing our MCMC and FIM results, we see there 
can be significant differences. Most parameters give results 
consistent to within an order (or two) of magnitude. The 
best agreement is for t p , which is largely uncorrelated with 
the other parameters. The widths for M., and l show 

more severe differences; these parameters show the tightest 
degeneracies. The two methods do show signs of slowly con- 
verging with increasing SNR, as expected. 

As a consistency check, to verify that the mismatch be- 
tween the FIM and MCMC results is a consequence of pa- 
rameter correlations, we calculated one-dimensional FIMs, 
only varying the MBH mass, and compared these to widths 
computed from MCMCs only sampling in mass. These were 
found to be in good agreement. The majority (~ 87%) have 
standard deviations consistent to within a factor of two; the 
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rest within an order of mag nitudej^S ome small difference is 
expected because of numerical error from calculating deriva- 
tives for the FIM by finite differencing. 



8.3 Scientific potential 

Having quantified the precision with which we could infer 
parameters from an EMRB waveform, we can now consider 
if it is possible to learn anything new. 

Of paramount interest are the MBH mass and spin. The 
current uncertainty in the mass is ctm. = 0.36 x 10 6 M© 
(~ 8%). There are few runs amongst our data set that are 



14 One differed by more than an order of magnitude, and also 
failed to fulfil the (one dimensional) MM criterion; this was a 
numerical problem in calculating the FIM. 



not better than this: it appears that orbits of a \i — 10 Mq 
CO with periapses r p < 13r g should be able to match our 
current observational constraints. However, the EMRB is an 
independent measurement, and so a measurement of compa- 
rable precision to the current bound can still be informative. 
Accuracy of 1% could be possible if r p < 8r g . 

The spin is less well constrained. To obtain an uncer- 
tainty for the magnitude of 0.1, comparable to that achieved 
in X-ray measurements of active galactic nuclei, it appears 
that the periapsis needs to be r p < H?V For smaller peri- 
apses, the uncertainty can be much less, indicating that an 
EMRB could be an excellent probe. The orientation angles 
for the spin axis may be constrained to better than 0.1 for 
r p < llr g . It may well be possible to learn both the direc- 
tion and the magnitude of the spin. This could illuminate 
the MBH's formation. 

We have no a priori knowledge about the CO or its or- 
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Figure 12. Distribution widths as functions of periapse r p and SNR p. Light blue is used for the standard deviation, red is the scaled 
50-percentile range and green is the scaled 95-percentile range: all three coincide for a normal distribution. Filled circles are used for 
converged runs, open circles for those yet to converge. The dotted line indicates the current uncertainty for M.; the dashed lines the 
standard deviation for an uninformative prior, and the solid lines the total prior range. 
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Figure 12 - continued Distribution widths as functions of periapse r p and SNR p. Light blue is used for the standard deviation, red 
is the scaled 50-percentile range and green is the scaled 95-percentile range: all three coincide for a normal distribution. Filled circles 
are used for converged runs, open circles for those yet to converge. The dotted line indicates the current uncertainty for M»; the dashed 
lines the standard deviation for an uninformative prior, and the solid lines the total prior range. 
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Figure 12 - continued Distribution widths as functions of periapse r p and SNR p. Light blue is used for the standard deviation, red 
is the scaled 50-percentile range and green is the scaled 95-percentile range: all three coincide for a normal distribution. Filled circles 
are used for converged runs, open circles for those yet to converge. The dotted line indicates the current uncertainty for M»; the dashed 
lines the standard deviation for an uninformative prior, and the solid lines the total prior range. 
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bit, so anything we learn would be new. However, this is not 
particularly useful information, unless we observe multiple 
bursts, and can start to build up statistics for the dynam- 
ics of the GC. Using current observations for the distance 
to the GC, which could be further improved by the mass 
measurement from the EMRB, it is possible to infer a value 
for the mass fi from £. This could inform us of the nature 
of the object (BH, NS or WD) and be a useful consistency 
check. A small value of £, indicating a massive CO, would 
be unambiguous evidence for the existence of a stellar mass 
black hole. 



<^ 1Q0 r 



9 EXTRA-GALACTIC SOURCES 

We have so far only been concerned with properties of bursts 
from our own galaxy. This is the best source for bursts be- 
cause of its proximity. A natural continuation is to consider 



EMRBs from other MBHs. Rubbo et al. (2006) suggested 



that LISA should be able to detect EMRBs originating from 
the Virgo cluster, altho ugh the detectable ra te may be only 



10 4 yr 1 per galaxy ( jHopman et al. 2007). Detect ability 



depends upon the mass of the MBH; higher masses corre- 
spond to lower frequency bursts, which are harder to detect. 

Checking our nearest neighbours, we find bursts from 
Andromeda (M31) would not be detectable. This is because 
of the large mass of the MBH M M3 i = (1.4±g;|) x 10 8 M o 
( Bender et aLl|2005) . However, its companion M32 is more 
promising. It has a li ghter MBH M M 32 = (2.5±0.5) x 10 6 M o 
( Ver olme et al.|2002 ). The trend between the periapse radius 
and SNR is shown in Fig. |14| The fit is again for orbits with 
/* = ^/ GMm32 Jrp <lxl0 _3 Hzto avoid the bucket of the 
noise curve. Bursts for a 1M© (10M©) can be detected with 
p > 10 if the periapse is smaller than 7r g (13r g ). 

The general behaviour is the same as for the GC. Bursts 
from the two MBHs can be compared using their character- 



Figure 14. Signal-to-noise ratio as a function of periapse radius 
for a pL = 1M© CO about the MBH of M32. The plotted points 
are the values obtained by averaging over each set of intrinsic 
parameters. The best fit line is log (p) = — 2.65 log(r p /r g ) + 3.05. 
This is fitted to orbits with r p > 18.8r g and has a reduced chi- 
squared value of % 2 jv = 1.26. 



istic frequencies /* and scaled SNR 



p*[h] = 



(- 



M 



R 



kpc J Vl0 6 M G 



-2/3 



P[h], 



(42) 



where R and M are the appropriate distances and masses 
for the two MBHs. These scalings can be determined from 
the quadrupole piece o f (|14| assuming a characteristic length 
scale set by r p . Figure [l5]shows the trend for both galaxies. 
The difference in sky position is largely washed out through 
the motion of the detector. 



M31 and M32 are at a distance of 770 kpc (Karachent- 



sev et al. 2004). It therefore seems unlikely that bursts 



could be observed from the Virgo cluster at a distance of 
#virgo « 16.5 Mpc ( |Mei et al.||2007 ). 

Triangulum (M33) is believed not to have an MBH. 
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Figure 13. Parameter standard deviations versus periapsis r p , showing dependence (or lack thereof) upon the spin magnitude |a*|. 



Merritt et al. (2001) use dynamical constraints to place 
an upper b ound on the mass of a central BH of Mm33 < 
3 x 1O 3 M , [Gephardt et alT] ( |200l| ) find a bound of M M3 3 < 
1.5 x 1O 3 M0. Observations of the ultra- luminous nuclear X- 
ray source (ULX) closest to the centre of M33 yield a best 
estimate of Mulx ~ O (10) Mq for the source object's mass 
( Foschini et~aLl|2004| [Weng et al.||2009| ). This is consistent 
with there being no MBH; the ULX originates from a stel- 
lar mass BH that is coincidentally located close to the core 



of the galaxy. Consequently, we do not expect to see any 
bursts from M33: to detect one would confirm the existence 
of a previously invisible MBH. 



10 DISCUSSION 

We have outlined an approximate method of generating 
gravitational waveforms for EMRBs originating at the GC. 
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Figure 15. Scaled signal-to-noise ratio as a function of charac- 
teristic frequency. 



This assumes that the orbits are parabolic and employs a nu- 
merical kludge approximation. The two coordinate schemes 
for a NK presented here yield almost indistinguishable re- 
sults. We conclude that either is a valid choice for this pur- 
pose. There may be differences when the spin is large and the 
periapse is small: ~ 10% for r p ~ 4r g , ~ 20% for r p c± 2r g . 

The waveforms created appear to be consistent with 
results obtained using Peters and Mathews waveforms for 
large periapses, indicating that they have the correct weak- 
field form. The NK approach should be superior to that of 
Peters and Mathews in the strong-field regime as it uses the 
exact geodesies of the Kerr spacetime. Comparisons with 
energy fluxes from black hole perturbation theory indicate 
that typical waveform accuracy may be of order 5%, but this 
is worse for orbits with small periapses and may be ~ 20%. 
These errors are greater than the differences resulting from 
the use of the alternative coordinate systems. 

The signal-to- noise ratio of bursts is well correlated with 
the periapsis. Except for the closest orbits (r p < 7r g ), the 
SNR (per unit mass) may be reasonably described as having 
a power-law dependence of 



log (p) ~ -2.7 log 



+ 4.9. 



(43) 



Signals should be detectable for a 1M© (10M©) object if 
the periapse is r p < 27r g (r p < 65r g ), corresponding to a 
physical scale of 1.7 x 10 11 m (4.1 x 10 11 m) or 5.6 x 10~ 6 pc 
(1.3 x 10~ 5 pc). 

Using the NK waveforms we conducted an investigation, 
using Fisher matrix analysis, into how precisely we could 
infer parameters of the galactic centre's MBH should such an 
EMRB be observed. However, we found that the linearised- 
signal approximation does not hold for these bursts over 
a wide range of SNR. This demonstrates the necessity of 
checking the approximation before quoting the results of an 
analysis (Vallisneri 2008). 

We used MCMC results as a more robust measure of 
parameter estimation accuracy. Potentially, it is possible to 
determine very precisely the key parameters denning the 
MBH's mass and spin, if the orbit gets close enough to the 
MBH. It appears that we can achieve good results from a 
single EMRB with periapsis of r p c± 10r g for a 1OM CO. 



This translates to a distance of 6 x 10 10 m or 2 x 10 -6 pc. 
Orbits closer than this would place stricter constraints. The 
best orbits yield uncertainties of almost one part in 10 5 for 
the MBH mass and spin, far exceeding existing techniques. 
Conversely, orbits with r p > 20r g are unlikely to provide 
any useful information. 

Before we can quote results for how accurately we can 
determine the various parameters, we must consider the 
probability of each orbit. This will be the subject of a com- 
panion paper, building upon the earlier results of Rubbo 
|et al.| ( |2006| and |Hopman et aL (2007), who only considered 
approximate forms for the SNR, rather than using wave- 
forms. Using a model for the nuclear star cluster of the GC 
it is possible to define distributions for angular momenta 
Loo, for a species of mass p. With these we can estimate the 
event rate and how much information, on average, we could 
hope to obtain from EMRB observations. If it is likely that 
we would observe multiple EMRBs, it may be possible to 
combine results to tighten uncertainties. 

Some consideration should also be given to meth- 
ods of fitting a waveform to an observed signal. Given a 
noisy data stream, how could EMRBs be extracted? The 
parabolic spectrum has a characteristic profile, suggest- 
ing that matched filtering could be possible. Complications 
could arise in fitting parameters to a waveform: we have 
seen that there exist complicated degeneracies between pa- 
rameters. These issues would warrant further investigation 
should the event rate be high enough. 

While we have only considered bursts from our own 
galaxy in detail, it should be possible to observe bursts from 
other nearby galaxies if their MBH is of the appropriate 
mass. This makes M32 a viable candidate. The SNR shows 
a similar dependence upon periapsis as for the GC, and may 
be described by a power-law of 



log (p) ~ -2.7 log 



+ 3.1, 



(44) 



for orbits with r p > 10r g . For a 1M© (10M©) object, bursts 
should be detectable for periapses r p < 7r g (r p < 14r g ), 
corresponding to 2.6 x 10 10 m (4.9 x 10 10 m) or 8.4 x 10~ 7 pc 
(1.6 x 10 -6 pc). This is a small region of parameter space, so 
we conclude that extra-galactic bursts are likely to be rare. 
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APPENDIX A: WINDOW FUNCTIONS 

When performing a Fourier transform using a computer we 
must necessarily only transform a finite time-span r. The ef- 
fect of this is the same as transforming the true, infinite sig- 
nal multiplied by a unit top-hat function of width r. Trans- 
forming yields the true waveform convolved with a sine. If 
h! (/) is the computed Fourier transform then 



h'U) 



Pt/2 

J h(t)e 27Vlft dt = \h(f) * r sinc(7r/r)J , (Al) 



where h(f) = & {h(t)} is the unwindowed Fourier transform 
of the infinite signal. This windowing of the data is a problem 
innate in the method, and results in spectral leakage. 



Figure 1(a) shows the computed Fourier transform for 
an example EMRB. The waveform has two distinct regions: 
a low- frequency curve, and a high-frequency tail. The low- 
frequency signal is the spectrum we are interested in; the 
high-frequency components are a combination of spectral 
leakage and numerical noise. The 0(1/ f) behaviour of the 
sine gives the shape of the tail. This has possibly been 



misidentified in figure 8 of Burko Sz Khanna] (12007} as the 
characteristic strain for parabolic encounters. 

Despite being many orders of magnitude below the peak 
level, the high-frequency tail is still well above the noise 
curve for a wide range of frequencies. It therefore contributes 
to the evaluation of any inner products, and could mask 
interesting features. It is possible to reduce the leakage using 
apodization: to improve the frequency response of a finite 
time series one can use a weighting window function w(t) 
which modifies the impulse response in a prescribed way. 

The simplest window function is the rectangular (or 
Dirichlet) window wn(t); this is the top-hat described above. 
Other window functions are generally tapered p*] There is a 
wide range of window functions described in the literature 



flHarris||1978l |Kaiser fc Schafer|[l980l |Nuttall|[l98l| |McK 
echan, Robinson, &; Sathyaprakash|2010 ). The introduction 
of a window function influences the spectrum in a manner 
dependent upon its precise shape. There are two distinct dis- 
tortions: local smearing due to the finite width of the centre 



15 When using a tapered window function it is important to en- 
sure that the window is centred upon the signal; otherwise the 
calculated transform has a reduced amplitude. 
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(a) Spectrum using no window. The calculated SNR is p = (b) Spectrum using a Nuttall window. The calculated SNR 
12.5. is p = 8.5. 



Figure Al. Example spectra calculated using (a) a rectangular window and (b) Nuttall's four-term window with continuous first 
derivative ( |Nuttall|198l| >. The spin of the MBH is a* = 0.5, the mass of the orbiting CO is p = IOMq, the periapsis is r p = 50r g and 
the inclination is t — 0.1. The high-frequency tail is the result of spectral leakage. The level of the LISA noise curve is indicated by the 
dashed line. The spectra are from detector I, but the detector II spectra look similar. 



lobe, and distant leakage due to finite amplitude sidelobes. 
The window function may be optimised such that the peak 
sidelobe has a small amplitude, or such that the sidelobes 
decay away rapidly with frequency. Choosing a window func- 
tion is a trade-off between these various properties, and de- 
pends upon the particular application. 

For use with the parabolic spectra, the primary concern 
is to suppress the sidelobes. Many windows with good side- 
lobe behaviour exist; we consider three: the Blackman-Harris 
minimum four- term window ( Harris|1978 Nuttall| 1981 ) 



n=0 



BH 



2n7Tt\ 



where 



BH 

a 2 



0.35875, 
0.14128, 



a? H = 0.48829, 



BH 

«3 



: 0.01168; 



(A2) 



(A3) 



the Nuttall four-term window with continuous first deriva- 
tive flNuttall|198T ) 

wn(£) = cos f^Hl \ ? (A4) 

n=0 ^ ' 



where 

= 0.355768, a? = 0.487396, 
a£ = 0.144232, a* = 0.012604, 



(A5) 



and the Kaiser-Bessel window ( [Harris 1978 Kaiser & 
Schafer|1980) 



wkb(£; /3) 



h [/VI - (2i/r) 2 ] 



HP) 



(A6) 



where l v {z) is the modified Bessel function of the first kind, 
and /3 is an adjustable parameter. Increasing f3 reduces the 
peak sidelobe, but also widens the central lobe. 

The Kaiser-Bessel window has the smallest peak side- 
lobe, but the worst decay (1//); the Nuttall window has 



the best asymptotic behaviour (l// 3 ); the Blackman-Harris 
window has a peak sidelobe similar to the Nuttall window, 
and decays asymptotically as fast (slow) as the Kaiser-Bessel 
window, but has the advantage of having suppressed side- 
lobes next to the central lobe. 

Another window has been recently suggested for use 
with gravitational waveforms: the Planck-taper window 



(Damour, Iyer, & Sathyaprakash 2000 McKechan et al 



2010) 



wp(t] e) 



with 
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exp(Z+) 
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+ 1 



^ t < 



exp(Z_) + 1 
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<t <r 



2e 



1 



(A7) 



(A8) 



_l±2(t/r) 1 -2e±2(t/r) 

This was put forward for use with binary coalescences, and 
has superb asymptotic decay. However, the peak sidelobe is 
high, which is disadvantageous here. We therefore propose 
a new window function: the Planck-Bessel window which 
combines the Kaiser-Bessel and Planck-taper windows to 
produce a window which inherits the best features of both, 
albeit in a diluted form, 



wp B (t]l3,e) = wp(t;e)w K B(t; f3). 



(A9) 



The window functions' frequency responses are plotted in 
Fig. |A2| There is no window that performs best everywhere. 

Figure |A1| shows the computed Fourier transforms for 
an example EMRB using no window (alternatively a rectan- 
gular window), and the Nuttall window^] Using the Nuttall 

16 The Blackman-Harris, Kaiser-Bessel and Planck-Bessel win- 
dows give almost identical results. 
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Figure A2. Window function frequency response. To avoid clut- 
ter, the response function is only plotted in detail until fr = 8, 
above this a smoothed value is used, as indicated by the dot- 
dashed line. As well as having good asymptotic behaviour, the 
Planck-taper window has the narrowest main lobe, except for the 
rectangular window. 



window, the spectral leakage is greatly reduced; the peak 
sidelobe is lower, and the tail decays away as 1/ f 3 instead 
of 1/ f. The low frequency signal is not appreciably changed. 

The choice of window function influences the results as 
it changes the form of h(f). The variation in results be- 
tween windows depends upon the signal: variation is great- 
est for low frequency bursts, as then there is greatest scope 
for leakage into the detector band; variation is least signifi- 
cant for zoom- whirl orbits as then there are strong signals to 
relatively high frequencies, and spectral leakage is confined 
mostly to below the noise level. To quantify the influence of 
window functions, we studied the diagonal elements of the 
Fisher matrix from a selection of orbits with periapses rang- 
ing from ~ 10r g -300r g . For orbits with small periapses all 
five windows (excluding the rectangular window) produced 
very similar results: the Planck-taper window differed by 
a maximum of ~ 0.5% from the others, which all agreed 
to better than 0.1%. The worst case results came from the 
lowest frequency orbits, then the Planck-taper window de- 
viated by a maximum of ~ 30% in the value for the Fisher 
matrix elements, the Blackman-Harris deviated by ~ 20% 
and the others agreed to better than ~ 5%. The Planck- 
taper window's performance is limited by its poor sidelobe 
behaviour; the Blackman-Harris has the worst performance 
at high frequency. 

For this work we have used the Nuttall window. Its 
performance is comparable to the Kaiser-Bessel and Planck- 
Bessel windows, but it is computationally less expensive as 
it does not contain Bessel functions. Results should be accu- 
rate to a few percent at worst, and results from closer orbits, 
which provide better constraints, should be less affected by 
the choice of window function. Therefore, we are confident 
that none of our conclusions are sensitive to the particular 
windowing method implemented. 

This paper has been typeset from a T^X/ IOT^X file prepared 
by the author. 



© 2012 RAS, MNRAS 000,|Tp5l 



